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Abstract Nuclear density functional theory (DFT) is one of the main theoretical tools used to study the 
properties of heavy and superheavy elements, or to describe the structure of nuclei far from stability. While 
on-going efforts seek to better root nuclear DFT in the theory of nuclear forces [see Duguet et al, this 
issue], energy functionals remain semi-phenomenological constructions that depend on a set of parameters 
adjusted to experimental data in finite nuclei. In this paper, we review recent efforts to quantify the related 
uncertainties, and propagate them to model predictions. In particular, we cover the topics of parameter 
estimation for inverse problems, statistical analysis of model uncertainties and Bayesian inference methods. 
Illustrative examples are taken from the literature. 

PACS. 21.60.Jz Nuclear Density Functional Theory and extensions - 21.10.-k Properties of nuclei; nuclear 
energy levels - 02.30.Zz Inverse problems 02.60.Pn Numerical optimization - 02.70.Uu Applications of 
Monte Carlo methods 


1 Introduction 

Applications of nuclear science in energy production or 
national security are based on nuclear data such as cross- 
sections, energy levels, and lifetimes. In many cases of in¬ 
terest, experimental measurements are not available, and 
guidance from theory is indispensable. In the valley of sta¬ 
bility, one still can employ simple models heavily tuned to 
existing data. For example, the fission model implemented 
in the GEF code uses about 50 parameters such as fission 
barrier heights and level density parameter which are pa¬ 
rameterized as a function of Z , N, or neutron incident 
energy. The code also uses databases of binding energies 
and shell corrections (in the ground state only). Based on 
these parameters and data banks, qualitative arguments, 
and a Monte Carlo sampling scheme, observables such as 
fission probabilities, fission fragment yields, and neutron 
multiplicities can be reproduced accurately in the actinide 
region (with a few exceptions) pQ. While such an empiri¬ 
cal approach fulfills some of the needs of data evaluators, 
however, its predictive power beyond the region where the 
model is fitted is null. Indeed, such models do not contain 
any physics principle related to nucleons, their interaction, 
and the quantum nature of the atomic nucleus. 

Therefore, even if data-driven empirical models will al¬ 
ways be helpful in the short term, one must try to root 
data evaluation into more microscopic theories of nuclear 
structure and reactions in order to gain confidence in the 


reliability of evaluations. In heavy elements, density func¬ 
tional theory is currently the only candidate for such a 
microscopic approach to nuclear structure. In particular, 
recent advances in high-performance computing have en¬ 
abled large-scale calculations of nuclear properties at the 
scale of the mass table |21I51 ■ Despite this progress, how¬ 
ever, the accuracy and precision needed in data evalua¬ 
tions represent a formidable challenge for nuclear density 
functional theory (DFT). As an example, nuclear bind¬ 
ing energies are computed within approximately 500 keV 
in state-of-the-art DFT calculations [115] . Although this 
represent a relative error of 0.05% or less for nuclei with 
mass A > 100, it remains far from the sub-keV accuracy 
that is demanded in, for example, criticality studies. In 
order to make data evaluations based on microscopic in¬ 
puts from DFT a viable alternative to simpler models, two 
challenges must be addressed in the next few years. 

First, DFT must be more firmly and rigorously con¬ 
nected to the theory of nuclear forces as defined, for exam¬ 
ple, by effective field theory [B] ■ One possibility is to derive 
local energy functionals from chiral effective field poten¬ 
tials by using the density matrix expansion 171151 . Another 
is to use, for example, many-body perturbation theory to 
expand energy and norm kernels of ab initio approaches 
in a form amenable to DFT treatment [5]; see also Duguet 
et al. in this issue. 

Second, irrespective of its mathematical form and phys¬ 
ical origin, DFT kernels will always contain a phenome- 
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nological component in the sense that they depend on a 
small set of parameters that must be adjusted to data. In 
addition to making the theory usable, adjusting these pa¬ 
rameters will effectively provide an ad hoc mechanism to 
capture missing correlations. However, this optimization 
will induce an obvious dependence on the data and the 
optimization process itself, in addition to the pre-existing 
uncertainties related to the form of the functional and pos¬ 
sible truncation errors in the numerical implementation. 
As a result of these uncertainties, it is by no means guar¬ 
anteed that the theory will be capable of converging to 
the exact result. Therefore, rigorous methodology is es¬ 
sential in order to identify, quantify and propagate model 
uncertainties. 

In this paper, we review the progress made in this area 
over the past 10 years. In sect. [2j we recall the essential as¬ 
pects of nuclear density functional theory; in particular we 
discuss the distinction between the self-consistent mean- 
field theory and the energy density functional (EDF) ap¬ 
proach. In sect. 3j we look at DFT from a statistician’s 
point of view: What are the parameters of the model? How 
can they be determined? How can we quantify the statis¬ 
tical uncertainties? We summarize the various tools used 
to answer these questions. In sect. [4j we review the most 
recent attempts to propagate statistical uncertainties in 
model predictions using either covariance techniques or 
Bayesian inference. 


2 Nuclear Density Functional Theory 

Density functional theory is a general approach for solving 
the quantum many-body problem. Its most rigorous for¬ 
mulation is in electronic structure theory, where it is based 
on the existence theorem of Hohenberg and Kohn [TO]. It 
states that the energy of an interacting electron gas can 
be written as a functional of the one-body local density 
(of electrons), and the minimum of this functional gives 
the exact ground-state of the system. Shortly thereafter, 
this formal existence theorem was supplemented with the 
Kohn-Sham scheme, which allows one to determine the ac¬ 
tual density of electrons that minimizes the energy (if the 
functional itself is known) by solving equations analogous 
to Hartree equations El- Various extensions have been 
proposed to handle exchange energy exactly (the Kohn- 
Sham equations then are similar to Hartree-Fock equa¬ 
tions), excited states, systems at finite temperature, and 
to reproduce superfluid correlations (see, e.g., 

These extensions rely on reformulating the Kohn-Sham 
scheme with the full one-body density matrix (rather than 
the local density), density operators, a combination of one- 
and two-body densities, and so forth. Often other exten¬ 
sions account for relativistic effects [IT] . 

Implementations of DFT in nuclear physics are less 
straightforward, since the nuclear Hamiltonian is not known, 
in contrast with electronic structure theory. In addition, 
nuclei are self-bound, and correlation effects are much 
stronger than in electron systems [6]. Consequently, most 
nuclear energy functionals used so far have been in fact de¬ 
rived from the expectation value on the quasiparticle vac¬ 


uum of effective nuclear forces used in the self-consistent 
nuclear mean-field theory m • Therefore, they are formu¬ 
lated in terms of the intrinsic one-body nonlocal density 
matrix and nonlocal pairing tensor, which can break sym¬ 
metries of realistic nuclear forces such as translational or 
rotational invariance, parity, time-reversal invariance, and 
particle number. This spontaneous symmetry breaking is 
essential for introducing long-range correlations in the nu¬ 
clear wavefunction rmfi7| . Nuclear energy functionals are, 
therefore, substantially different from their counterpart in 
electronic DFT. This difference is reflected in the appel¬ 
lation of energy density functional (EDF) formalism. 

2.1 The Energy Density Functional Approach 

In this section, we succinctly describe the basic ingredients 
of the single-reference EDF (SR-EDF) approach with non- 
relativistic empirical functionals such as derived from the 
Skyrme or Gogny effective interactions. We refer to [61117] 
for discussions of more general frameworks such as mul¬ 
tireference EDF and ab initio DFT. The starting point is a 
set of single-particle states |*) that form a basis of the one- 
body Hilbert space. The related creation/annihilation op¬ 
erators are c\ and Cj and define the configuration space 
representation of the Fock space |18|fT6|. The coordinate 
space representation is obtained by invoking the continu¬ 
ous basis |ir) = | rar) of the one-body Hilbert space, with 
a = ±1/2 the intrinsic spin projection and r = ±1/2 the 
isospin projection. The single-particle functions are then 
(x\i) = (j>{x), and the corresponding creation/annihilation 
operators are the field operators c'(x) and c(x). The parti¬ 
cle vacuum of the Fock space is denoted by |0) and is char¬ 
acterized by Vi, Ci|0) = 0, or, alternatively, Vx, c(ir)|0) = 
0. 

Because of the importance of pairing correlations in 
low-energy nuclear structure m, we introduce a canoni¬ 
cal transformation between particle operators and quasi¬ 
particle operators /3^,/^. This Bogoliubov-Valatin trans¬ 
formation is mrnmm 

^ ' [^/im Gn ± k/m Gn] ’ 

m (1) 

^ = J2KmC m +U^ m cl] . 

m 

In the SR-EDF approach, we introduce the reference state 
| <P) as a product wavefunction of quasiparticle operators 
acting on the particle vacuum, 

i*>=n«. ( 2 ) 

Note that, by construction, the quasiparticle vacuum ([ 2 ]) 
does not conserve the particle number. 

The next step is to recall that for any given many-body 
state |<F), the one-body density matrix p and two-body 
pairing tensor n are defined in configuration space as 

_ <p\c\a\*) _ WcjCiW 

Pij (V\&) ’ Kij (V\V) 


(3) 
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and the generalized density 1Z as 


U = 



(4) 


When | E) is the quasiparticle vacuum Q, the correspond¬ 
ing generalized density matrix verifies 1Z 2 = 7Z and 1Z' = 
1Z jT8] . In addition, the Wick theorem ensures that the ex¬ 
pectation values of any operator on the quasiparticle vac¬ 
uum can be expressed as functions of p, k and k* alone. 
These three mathematical objects are thus the basic de¬ 
grees of freedom of the theory. In particular, the energy is 
then expressed as a functional E[p, k , «*]. 

The actual density and pairing tensor of the nucleus 
in its ground-state are determined by solving the Hartree- 
Fock-Bogoliubov (HFB) equations, which are obtained by 
applying the variational principle with respect to p, k, and 
k* [THUS]. This leads to 


[H,K]=0, (5) 

where 7~L is the HFB matrix, 7~L,j = dE/dlZji. One-body 
observables can then be computed as the trace of the rele¬ 
vant operator and p. Because of the nonlinear nature of the 
HFB equations, it is possible for the generalized density 
to break various symmetries of nuclear forces. Conversely, 
conserved symmetries can be used to label quasiparticle 
states msmm . 


2.2 Pseudopotentials and Energy Functionals 


Until recently, most applications of the nuclear EDF ap¬ 
proach have been based on semi-empirical EDFs explicitly 
derived from the expectation value of effective two-body 
forces U e ff. on the quasiparticle vacuum, 


E[p,k,k*] 


($\T + V e s.\$) 


( 6 ) 


In particular, the Skyrme effective force is a zero-range 
two-body pseudopotential for which the EDF becomes a 
functional of the local density only [251126] . The Gogny 
force has a finite range and gives a functional of the non¬ 
local one-body density EZ|; see, for example, [~P7f28l for 
comprehensive reviews of applications of Skyrme and Gogny 
EDFs. The empirical nature of both the Skyrme and Gogny 
potentials is manifested by the presence of density depen¬ 
dencies, which prohibits writing the potential in strict sec¬ 
ond quantization form [29] . With the exception of a few re¬ 
cent applications mmm , these EDFs have been used in 
the context of the self-consistent mean-field theory rather 
than in a strict Kohn-Sham scheme. 

In particular, many applications used the underlying 
effective pseudopotential V e e. to implement beyond mean- 
field techniques, where EDF reference states of the type 
([2| serve as basis states to expand the unknown rnany- 
body wavefunction, for example, in the generator coor¬ 
dinate method, or to restore broken symmetries by us¬ 
ing projection techniques mm A few years ago, how¬ 
ever, standard beyond mean-field techniques were shown 


to be invalid with density-dependent pseudopotentials [33l 
nmim This result has stimulated efforts to remove 
density dependencies, for example by using momentum- 
dependent two-body pseudopotentials [38] or zero-range 
two- and three-body pseudopotentials mm- Since these 
pseudopotentials are specifically designed to enable be¬ 
yond mean-field techniques such as projection and config¬ 
uration mixing, the central element of all these approaches 
is the effective Hamiltonian T + V e s. rather than the EDF 
itself. 

An alternative route is to implement a strict Kohn- 
Sham approach, where the only degrees of freedom are 
p, k, and k* , the ground-state wave function is always a 
quasiparticle vacuum of the form ([2]), and there is no men¬ 
tion of some underlying effective potential V e s, . In such 
an approach, the energy functional E[p, k, k*] must be de¬ 
signed so that it contains all relevant types of correlations. 
Only two main families of such functionals have been pro¬ 
posed in the literature: those proposed by Fayans and 
collaborators mmm, and the BPCM functionals from 
the Barcelona-Paris-Catana-Madrid collaboration FTTlFm . 
The main difficulty of this strict Kohn-Sham scheme, which 
is more in line with the spirit of DFT as encountered in 
electronic structure theory, is to incorporate beyond mean- 
field correlations accounting, for example, for large ampli¬ 
tude collective motion, or symmetry restoration. Recent 
work suggests that this could be achieved by introducing 
new densities representing collective degrees of freedom 
such as two-body or “collective” densities [ 17 ) 113711331133 ] 
(which may lead to a generalization of the Kohn-Sham 
equations) or by adding specific terms to the functional 
designed to cancel symmetry-breaking mm- 


2.3 Pairing Correlations 


Sometimes overlooked is the fact that, according to the 
HohenbergKohn theorem, the exact ground-state of the 
system can in principle be expressed entirely as a func¬ 
tional of the local one-body density matrix only. If this 
theorem could be extended directly to nuclear function¬ 
als of the intrinsic one-body density (see, e.g., mm, 
pairing correlations could — in principle — be produced 
by an unique functional of p(r). In this idealized scenario, 
there would be no need for quasiparticle operators, the Bo- 
goliubov transformation, or the pairing tensor: the Kohn- 
Sham scheme would be implemented directly with EDF 
reference states taken as particle number conserving Slater 
determinants. 

In practice, of course, the form of this functional is 
totally unknown. Until further notice, therefore, it seems 
more reasonable to build on the success of the self-consistent 
mean-field theory, to seek an explicit pairing term that is 
a functional of the usual pairing tensor, and to work with 
symmetry-breaking reference states of the type ([2]). Recall 
that the pairing tensor is defined from the specific form 
that the two-body correlation function takes in an HFB 
vacuum mm- If we denote p 2 {%i, X 2 , x[, x' 2 ) as the full, 
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nonlocal, two-body density matrix, then we have 

p 2 {x 1 ,X2,x' 1 ,x' 2 ) = K*( X 1 i X 2 )k(x' 2 ,X I 1 ) 

- p{x' 2 ,x 1 )p(x[,x 2 ) + p(x' 1 ,x 1 )p(x' 2 ,x 2 ). (7) 

Because of this property, the pairing tensor re and its com¬ 
plex conjugate re* are, indeed, the two only degrees of free¬ 
dom needed to account for pairing correlations at the HFB 
level. 

The pairing EDF can then be obtained by taking the 
expectation value of a pseudopotential V^ air ^ on the quasi¬ 
particle vacuum, which will immediately introduce a de¬ 
pendence on re* and re. This potential can be the same as 
the one used in the particle-hole channel, which is typi¬ 
cally the choice retained when working with the Gogny 
force m■ It can also have a different form, ranging from 
simple seniority pairing forces (16} to density-dependent 
zero-range pairing forces [54} to separable expansion of 
finite-range, Gogny-like potentials [551156] . Most of these 
pairing forces, and hence the resulting pairing functionals, 
are characterized by only a few parameters, and all lead 
to EDFs that are functionals of re* and re only. 


3 Density Functional Theory as a Model 


7]ti{x) is thus the output of an HFB calculation for the ith 
observable of type t. Also, e is the vector containing the 
error between the actual calculation and the experimental 
value. By definition, we thus have 


Uti = Vti{x) + e u , e ti ~ iV(0,cr t ), V(f,t). (8) 

Based on these notations, we minimize the weighted mean 
squared deviation given by 


X 2 {x) = 


n d - n, 


EE 

t =i i =i 


yu - vti{x) 




(9) 


where T is the total number of different data types, nt the 
total number of points of type f, and n d the total number 
of experimental points, n d = XEi n t ■ By convention, the 
vector of parameters at the minimum of the y 2 is noted 
x. Recall that y 2 1 implies a poor fit, where y 2 ss 1 
indicates a good fit. This is the familiar “y 2 per degree of 
freedom.” If all errors e t i are independent and normally 
distributed with mean 0, then the minimization of (|9| ) is 
equivalent to maximizing the likelihood function |59ll60j . 
In addition, the y 2 is a random variable that follows a 
genuine y 2 probability distribution function. 


Whether building the description of atomic nuclei on an 
effective potential U e g that defines both the mean-field 
and beyond mean-field corrections, or an EDF E\p, re, re*] 
in a strict Kohn-Sham framework, theoretical predictions 
will depend on a set of unknown parameters x correspond¬ 
ing, respectively, to the parameters of the effective nuclear 
force or the coupling constants of the EDF. Some of these 
parameters may be constrained by exploring the connec¬ 
tions with the theory of realistic nuclear forces or inves¬ 
tigating ideal systems such as nuclear matter or neutron 
drops [571155] . In general, however, one will also have to 
introduce experimental data in nuclei in order to set the 
values of these parameters. This fit of low-energy coupling 
constants to experimental data belongs to the class of in¬ 
verse problems in statistics. In this section, we review some 
of the techniques used in nuclear DFT to solve this prob¬ 
lem. Most of our considerations are based on the SR-EDF 
approach to nuclear structure but are easily extended to 
the self-consistent mean-field approach. 

3.1 Parameter Estimation 

The problem of determining the parameters of the nuclear 
EDF is easily posed: one needs only to choose a set of data 
points, define an objective function such as a y 2 function, 
and minimize the objective function with respect to the 
parameters. We use the following notations: y denotes the 
values of a set of experimental observables, with yu the 
value of the ith observable of type t; x = {x\ ,..., x na .) 
represent the vector of the n x parameters of the model, 
that is, the EDF coupling constants in our case; and r) 
collects the output of all model calculations. In our case, 


3.1.1 Experimental Dataset and Bias Estimation 

Choosing which and how many data points to include in 
the y 2 is the first important decision, and several strate¬ 
gies have been followed. In nuclear mass models based 
either on the Skyrme or Gogny force, all available experi¬ 
mental information on atomic masses is used; see [5114] and 
references therein. It is supplemented by additional data 
on, for example, fission barriers [fill or neutron matter [4]. 
The main concern for mass models is the risk of producing 
a high-bias estimator of the data. In simpler terms, it is 
by no means guaranteed that mass model parameteriza- 
tions of Skyrme of Gogny forces are reliable for computing 
observables that are not masses. 

By contrast, most historical fits of the Skyrme and 
Gogny forces were based on the smallest possible set of 
data. These included nuclear matter properties, binding 
energies, radii, and single-particle states; see [T511281I62] for 
a discussion. In addition, data in finite nuclei was taken 
almost exclusively in doubly-magic spherical nuclei. Two 
of the most notable exceptions are the SkM* parameteri¬ 
zation of the Skyrme force [55] and the D1S parameteriza¬ 
tion of the Gogny force |64j . which included information 
on the fission barrier in 240 Pu. The combination of small 
dataset and spherical nuclei is also likely to lead to high- 
bias estimators. 

The recently proposed UNEDF parameterizations of 
the Skyrme EDF represent an attempt to reduce the bias 
of the fitting procedure by selecting a medium-sized sam¬ 
ple of 100+ data points carefully selected in both spherical 
and deformed nuclei; see fig. [T| for the specific case of the 
UNEDF2 functional [551I5T1I52] . As a result, the ability of 
UNEDF functionals to reproduce masses or fission barriers 
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has degraded between UNEDFO (3 datatypes, 108 points) 
and UNEDF2 (5 data types, 130 points), while the pre¬ 
dictive power on binding energies near closed shells and 
on single-particle states increased. Note that until now, no 
attempt has been made to rigorously quantify the bias of 
EDF parameterizations. 
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Figure 1. Experimental dataset used for optimizing the UN- 
EDF2 Skyrme functional [321 . 

The pairing channel represents an additional difficulty 
when determining the parameters of the nuclear EDF. In¬ 
deed, very little data can effectively and unambiguously 
constrain the pairing functional directly at the HFB level. 
In practice, the odd-even staggering of binding energies is 
most often used [551[MlfB71l68i[M] . The UNEDF function¬ 
als were the first ones where the fit of the pairing func¬ 
tional was performed simultaneously with the fit of the 
Skyrme EDF. As a result, there are built-in correlations 
between the parameters of the Skyrme EDF and the two 
parameters that control the pairing functional. 

3.1.2 Optimization Algorithm 

The minimization of the x 2 function in the context of nu¬ 
clear DFT remains costly in computational resources. In 
the example of the UNEDF functionals, 100+ full, axially- 
deformed HFB calculations must be performed in order 
to define the \ 2 - Some of the most popular parameteri¬ 
zations of the nuclear EDF were published in the 1980s 
and 1990s, where the cost of running a full HFB calcu¬ 
lation was prohibitive in terms of y 2 minimization. Even 
now, the optimization of HFB mass models, where the x 2 
includes over 2,500 points, is still performed in spherical 
symmetry by using empirical renormalizations ma¬ 
in view of this computational cost, specifically designed 
algorithms with a focus on efficiency and robustness are 
especially valuable. We recall that derivatives drj t i(x)/dx M 
are not available analytically for the minimization of the 
X 2 Of course they can always be computed numeri¬ 
cally, but at a significant cost when n x is large. Therefore, 


the optimization of the nuclear EDF is most efficiently 
performed with derivative-free approaches. 



Number of Evaluations 


Figure 2. Performance (log-log scale) of three solvers (limited- 
memory variable metric, POUNDERS Nelder-Mead) for non¬ 
linear generalized \ 2 problems with n x = 6, rid = 428; from 

ED- 

Although a rich literature on the subject exists, we 
mention here only the POUNDERS algorithm developed in 
the framework of the UNEDF collaboration |3nU721i73U7T] . 
POUNDERS is a derivative-free trust-region method based 
on forming a local quadratic model of each component of 
the x 2 - The quadratic models are valid only in a small 
region of the parameter space near the current point x , 
but their aggregate approximation can be minimized ana¬ 
lytically. The minimum x + Sx defines the next point of a 
Newton-like procedure. As seen in fig. [2j the POUNDERS 
algorithm converges quickly compared to the traditional 
Nelder-Mead algorithm; most important, it gives a better 
solution. POUNDERS has recently been applied to other 
problems of interest in nuclear physics mm- 


3.1.3 Choice of the Objective Function 

Given a set of data points and a minimization algorithm, 
some latitude remains in defining the weights, or stan¬ 
dard deviations a t , associated with each data type. These 
quantities represent the estimated error on the data type 
t. They are in principle determined in such a way that 
the x 2 objective function § approaches 1 at the min¬ 
imum. Satisfying this condition may require readjusting 
the weights during the minimization HZ]. In practice, this 
step has rarely been done in nuclear EDF optimization, 
and the weights are most often kept constant (though data 
type dependent). Results reported in [76] and listed in ta- 
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Table 2. Rerun of POUNDERS on the UNEDFO problem (rid = 108, njv = 72) from two different starting points. The scaled 
difference columns represent the difference between the final value found and the original UNEDFO parameterization, scaled by 
its uncertainties cr,;; see table VIII in |30l . 



Starting from SLy4 

Scaled Diff. 

Starting from SKM* 

Scaled Diff. 


initial 

final 


initial 

final 


Pc 

0.159539 

0.160486 

-0.03954 

0.160319 

0.160435 

-0.09106 

E NM /A 

-15.9721 

-16.0685 

-0.2285 

-16 

-16.073 

-0.3119 

j^NM 

229.901 

230 

- 

216.658 

230 

- 

a NM 

^sym 

32.0043 

31.3393 

0.2604 

30.0324 

31.7221 

0.3856 

r NM 
-t^sym 

45.9618 

54.2493 

0.2290 

45.7704 

60.4725 

0.3844 

i/, m; 

1.43955 

0.9 

- 

1.26826 

0.9 

- 

,-ipAp 

°0 

-76.9962 

-55.2344 

0.01545 

-68.2031 

-55.7348 

-0.2794 

QP&P 

15.6571 

-64.1619 

-0.1499 

17.1094 

-70.4274 

-0.2599 

V 0 n 

-285.84 

-170.796 

-0.2003 

-280 

-170.788 

-0.1966 

v 0 p 

-285.84 

-197.782 

0.4238 

-280 

-198.038 

0.3474 

nP V J 
°0 

-92.25 

-77.9436 

0.4637 

-97.5 

-79.2915 

0.06990 

£rpV J 

-30.75 

27.4519 

-0.6171 

-32.5 

49.5737 

0.1339 

m 

1188.75 

67.9034 


24814.1 

67.5738 


n f 


235 



150 



Table 1. Root-mean-square deviations for each observable in 
the UNEDF1 optimization protocol compared for UNEDF1- 
HFB for a few different values of the standard deviation ao es 
(in MeV) of the OES data. All r.m.s. values are in MeV except 
the ones for proton radii, which are in fm; from Eg. 


CTOES 

0.025 

0.050 

0.075 

0.100 

Deformed masses 

0.944 

0.776 

2.596 

0.806 

Spherical masses 

2.427 

1.836 

2.669 

1.718 

Proton radii 

0.022 

0.022 

0.022 

0.022 

OES neutrons 

0.012 

0.051 

0.065 

0.080 

OES protons 

0.043 

0.074 

0.075 

0.072 

Fission isomer 

0.809 

0.558 

0.535 

0.530 


ble[T]show that the impact of the weights on the result of 
the optimization could be significant. 

Minimizing the y 2 § requires initializing the opti¬ 
mization algorithm with a vector Xq- Ideally, the optimiza¬ 
tion algorithm would be able to converge to the absolute 
minimum of the objective function, given a set of con¬ 
straints dictated by reality. In practice, it is nearly impos¬ 
sible to guarantee such a result. To our knowledge, there is 
only one example where the impact of the initial point on 
the resulting parameterization was studied in detail m- 
Table [2] illustrates the robustness of the POUNDERS algo¬ 
rithm; similar solutions to the optimization problem are 
obtained when starting from the SLy4 or SkM* parame¬ 
terization. The largest difference occurs for L ™f (slope of 
the symmetry energy in nuclear matter at saturation den¬ 
sity) and C’f VJ (isovector spin-orbit coupling constant), 
which are poorly constrained by the dataset. 


3.2 Statistical Uncertainties of Energy Densities 

Irrespective of the form of the energy functional, the de¬ 
gree of arbitrariness in defining the y 2 used to determine 


the “best” parameters of the EDF clearly suggests possi¬ 
bly large uncertainties in the resulting parameterizations. 
Standard methods of probability and statistics can be used 
to quantify some of these uncertainties. In this section, we 
review only the techniques used to estimate statistical un¬ 
certainties. Few studies of systematic uncertainties have 
been conducted so far; see m for discussion. Numerical 
errors are discussed separately in sect. |3.3| This separa¬ 
tion is made for convenience only, since it is illusory to 
think that all sources of uncertainties can be completely 
disentangled. 

3.2.1 Covariance Analysis 

One of the most common quantities used for estimating 
statistical uncertainties is the covariance matrix. The use 
of covariance techniques in nuclear DFT is relatively re¬ 
cent. Full regression analysis was first introduced in the 
context of nuclear mass fits in m The covariance ma¬ 
trix was first mentioned and computed for Skyrme EDF 
optimization in [62] . Since then, there have been many 
applications of this technique to compute the standard 
deviation of EDF parameters and propagate uncertainties 
in model predictions; see sect. [4] 

In the following, we denote Cm the covariance matrix 
of the parameters x of the model M that we are using 
(EDF), formally, 

( C M )ij = E [{Xi - E(xi))(xj - E(xj ))\, (10) 

where E() refers to the average of a random variable; each 
parameter Xi is thus treated as a random variable. One 
should distinguish Cm from the “data” covariance matrix 
Cd- The latter notation will be used to refer to the covari¬ 
ance matrix of the random variables e associated with the 
error between the model output Tj(x) and the experimen¬ 
tal data y. All these misfits eu are often assumed to be in¬ 
dependent, therefore Co is diagonal and (Co)ij = crfdij. 
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One also assumes that they follow a (multivariate) normal 
distribution with mean 0, e ~ AF(0, Cp). 

In the simple case of an unweighted, linear least-squares 
optimization, where rj{x) = Ax and of = 1, the inverse 
of the covariance matrix Cm can be computed as mm 


random variables, it reads 


p{x\y, M)dx 


P(y\x, M)P(x\M)dx 
[ P(y\M)dx 


(13) 


= dn-%, ( 11 ) 

where Ti is the Hessian matrix of (rid — n x )y 2 (x). In the 
case of nuclear EDF optimization, the quality of the co- 
variance matrix estimation is thus contingent on the linear 
dependence of observables with model parameters within 
the range of variation of interest. From the literature, one 
finds that nuclear binding energies behave linearly across 
a broad range of parameter space eg; single-particle or¬ 
bitals have a small degree of nonlinearity El; nonlinearity 
becomes more pronounced in the variation of fission iso¬ 
mer excitation energies m- While covariance techniques 
have been often employed recently to obtain estimates of 
statistical uncertainties on model predictions, see sect. [4j 
the underlying hypothesis of linearity has rarely been in¬ 
vestigated in detail. 

The covariance matrix can also be used to get an es¬ 
timate of confidence intervals/regions. Recall that if the 
errors e follow a multivariate normal distribution, then the 
confidence interval at a percent for parameter i is defined 
by the endpoints 

-G i \j [C Al putrid— iix, 1— ^ ; (12) 

where f nd _ nxi i_| is the 1 — f quantile of the (Student’s) 
t distribution with rid — n x degrees of freedom [8UJEM1E1 
[811182] , This was used, for example, in the assessment of 
the UNEDF functionals j30U3Hf5?| . The diagonal elements 
of the covariance matrix define the standard deviations, 
(Cm)u, of each parameter i. 

3.2.2 Bayesian Techniques 

Bayesian inference techniques have been used for many 
years in the nuclear data community pilMlf^llMI^T] . In 
nuclear structure, this method has recently gained ground, 
for example, to quantify uncertainties in chiral effective 
potentials (8811891 . In DFT, Bayesian inference has been 
used in electronic structure theory to evaluate uncertain¬ 
ties induced by the fit of the exchange-correlation func¬ 
tional in the generalized gradient approximation (90]. 

Following [83], one may describe Bayesian techniques 
as an exercise of inductive inference when the probability 
for an hypothesis A to be true is interpreted not strictly as 
the number of observations of A over the total number of 
outcomes, but rather as the degree of plausibility that A 
is true. The “philosophical” interpretation is that the true 
value of the model parameters x is described probabilis¬ 
tically. Additional data further constrain the probability 
distribution but never reduces it to a single, known value. 

The Bayes theorem provides the mathematical foun¬ 
dation for Bayesian techniques. In the case of continuous 


In practice, we seek the probability of the model M hav¬ 
ing parameters x based on a set of observed data y. The 
model is typically characterized by a number of features 
that add to the resulting uncertainties. In the case of EDF 
optimization, these features include the type of functional 
(Skyrme, Gogny, or other), the treatment of pairing corre¬ 
lations (HFB approximation, particle number projection), 
and the numerical implementation. The probability distri¬ 
bution p(x\y, M) is the posterior distribution. Note that 
the posterior is computed within the model M: it does 
not contain any information about the validity of said 
model. In other words, suppose the posterior distribution 
is sharply peaked around a given value x 0 : the fact that 
Xq is the most likely parameter set does not mean (i) that 
it is the correct one (since more data may change the dis¬ 
tribution), and (ii) that the resulting model is the correct 
one (since everything is model-dependent). 

In eq. P{y\x 1 M) is the probability that the model 
produces the data given the parameters: it is the likelihood 
function [59i. P(x\M) is the probability that the model 
has parameters x irrespective of any data: it is the prior 
distribution. For a uniform prior distribution, maximiz¬ 
ing the posterior distribution is equivalent to maximizing 
the likelihood. We remark in passing that both statistical 
approaches give different results if one looks at the prob¬ 
ability distribution of some new parameter g(x) that is a 
function of the original parameters x [ 591 . 

Bayesian inference can also be used to compute an 
estimate of the covariance matrix Cm between the pa¬ 
rameters. Assuming weak nonlinearities of the model pa¬ 
rameters, that is, ri(x) oc x, the likelihood function is ap¬ 
proximately Gaussian with respect to x. If one assumes, 
for simplicity, full ignorance about the prior distribution 
(uniform distribution with independent parameters), then 
the posterior covariance matrix is given by (91] 

C M 1 « G t C~ d 1 G , (14) 

where 

G ij( x ) = ^r( x ) ( 15 ) 

and Cd is, as before, the covariance matrix associated 
with the misfits between data and the predictions. Owing 
to the (near) linearity of the model parameters, one can 

easily find that C M — 2“H as obtained from the stan¬ 
dard covariance matrix. The advantage of the Bayesian 
approach is the possiblity of including in the calculation 
of the covariance matrix the effect of prior knowledge of 
the distribution of model parameters; see sect. 3.2.3 in [9l] 
for details. 

Posterior distributions are typically sampled by using 
Markov chain Monte Carlo (MCMC) techniques [S2], the 
result being a (dependent) sequence of samples {x^,..., x (2 )}. 
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In practice, this sampling can be computationally chal¬ 
lenging, since thousands or millions of evaluations of the 
likelihood function, he nce of the x 2 function, may be needed. 
As mentioned in sect. 3.1.2 the x 2 functions used in nu¬ 


clear EDF optimizations may typically involve between 
100 and 2,500 HFB calculations or more, making the di¬ 
rect sampling of the posterior distribution prohibitive. The 
alternative is to estimate response surface functions in or¬ 
der to emulate the behavior of the model response r]u(x) 
at a much cheaper cost mmm- The parameters of these 
response functions can also be incorporated in the statis¬ 
tical setting x. 
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only on the problem of solving the HFB equations: in the 
SR-EDF approach, these are the only equations needed. In 
multireference EDF, the HFB equations also play a cen¬ 
tral role because errors in the solutions will propagate to 
the calculation of beyond mean-field corrections such as 
in the GCM ii l5) . 


Mesh Size [fm] 



Figure 3. (Color online) Univariate and bivariate marginal 
estimates of the posterior distribution for the 12-dimensional 
DFT parameter vector of the UNEDF1 parameterization. The 
blue lines enclose an estimated 95% region for the posterior 
distribution found when only the deformed masses from the 
original UNEDF1 data are accounted for. The red dot corre¬ 
sponds to the UNEDF1 values; see [76] . 


Until now, there have been only two examples of Baye¬ 
sian applications in nuclear EDF optimization. In 196] , the 
backward forward Monte-Carlo method was applied to es¬ 
timate uncertainties in Skyrme mass model parameters. 
In [SZ], the full posterior distribution of the Skyrme EDF 
corresponding to the UNEDF1 x 2 was computed by using 
response functions based on Gaussian processes. The re¬ 
sulting 12-dimensional multivariate distribution is shown 
in fig. [3] The characteristics of the posterior distribution, 
such as the calculated standard deviations of the parame¬ 
ters, are similar to the results from the analysis based on 
confidence interval given in m- 


3.3 Numerical Implementations 

Implementing the DFT equation (e.g., the HFB equations 
for the SR-EDF approach) in a computer program intro¬ 
duces numerical errors. These errors are unavoidable be- 
caues the density matrix and pairing tensor have an infi¬ 
nite number of degrees of freedom. In this section, we focus 


Figure 4. Comparison between the pace of convergence of a 
spherical DFT calculation in coordinate-space, (red squares), 
and configuration space (HO basis), (black circles). Results 
were obtained by setting both direct and exchange terms of 
the Coulomb potentials to 0. The HO basis results are op¬ 
timized with respect to the oscillator frequency. Coordinate 
space calculations were performed with HFBRAD in a box of 
20 fm [98], HO calculations with HOSPHE [8]; from 1761 . 


One of the most popular approaches to solve the HFB 
equations is to expand the HFB solutions on a basis of 
known functions. In atomic nuclei, the eigenstates of the 
harmonic oscillator (HO) are most often used, since they 
are given analytically on spherical, cylindrical, and Carte¬ 
sian coordinates. In addition, there is an exact separa¬ 
tion between center of mass and relative motion. The 
lower part of the nuclear mean-field is well approximated 
by an HO. Several DFT solvers using HO basis expan¬ 
sions have been published; see [99lll00li22lll01l[l02]ll03( l50 1 
milling. In practice, all basis expansions are truncated. 
Therefore, HFB solutions become dependent on the char¬ 
acteristic parameters of the basis; in the case of the HO, 
these are the basis frequencies u> = (u> x ,u> y ,uj z ), number of 
oscillator shells N, and total number of basis states (if the 
basis is spherical, all frequencies are identical and there is 
a one-to-one correspondence between number of shells and 
number of states). This spurious dependence may induce 
large errors, for example in nuclei with large elongations 
or weakly bound systems. PMEHI- 

The HFB equations can also be solved by direct numer¬ 
ical integration; see, for example, una for the HFB for¬ 
malism in coordinate space. This has been done in spher¬ 
ical and axial symmetry only [9811108] . For more complex 
geometries, the computational cost of direct integration 
becomes prohibitive; and hybrid strategies such as lattice 
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discretization [10911110] . finite element analysis ITTT1HT21 
and multi-resolution wavelet expansion m have been 
investigated. While numerically significantly more precise 
(see fig. [I] for a comparison between the pace of conver¬ 
gence of the two methods), r-space-based techniques come 
with a higher computational cost, in terms of processes, 
memory, or disk space. Such approaches are also not ideal 
for handling finite-range local forces or nonlocal forces. 

Numerical errors inherent in DFT solvers are often 
overlooked, even though they may play a non-negligible 
role in the estimation of statistical uncertainties. For ex¬ 
ample, the truncation error of HO expansions increases 
with nuclear deformation, even when one tries to adjust 
accordingly the geometry of the HO basis [11411106] . As a 
result, the numerical error in the energy of, say, the fis¬ 
sion isomer or the top of the fission barriers in actinide 
nuclei is always going to be larger than the error in the 
ground state. In fact, at very large deformations, the er¬ 
ror of one-center basis expansions can reach a few MeV. 
Apart from adopting empirical corrections based on aux¬ 
iliary large-scale surveys of numerical errors m, the so¬ 
lution could be to generalize asymptotic formulas such as 
proposed in the context of ab initio theory |1161lll7llll8j . 
This problem, as well as the inclusion of these errors in 
the calculation of uncertainties, remains open. 


uation implies that computed observables are linearly de¬ 
pendent on model parameters, which is guaranteed only 
locally near the optimal point. The computed value r/ y (x) 
of a single new observable y depends on the parameter¬ 
ization of the EDF, and one can estimate its standard 
deviation based on the parameter covariance matrix Cm 

mM\- 


= Y. G y^ c MhG VJ 


G 


dr lv, 


yi(x) — ( X ), 


(16) 


If one now considers two new observables y and y ', pos¬ 
sibly correlated, such as the neutron skin in 208 Pb and 
electric dipole (El) polarizability an in the same nucleus, 
then the above formula should be generalized to 


Cw = G 1 C m 'G 


(17) 


to account for cross-correlations. 

In the context of DFT applications, such covariance 
analysis has been applied to compare statistical and sys¬ 
tematic uncertainties of neutron skins 11191 ; to explore the 
properties of ground-state properties of closed-shell nuclei 
far from stability and to optimize EDF for nuclear 

astrophysics [121(11221.1128] . 


4 Uncertainty Propagation and Predictive 
Power 

One of the main advantages of using the statistical anal¬ 
ysis techniques briefly presented in sect. [3] is to provide a 
rigorous framework for propagating the quantified uncer¬ 
tainties to predictions. These predictions can be the result 
of running the same model on a different dataset; for ex¬ 
ample, computing masses of exotic neutron-rich nuclei or 
superheavy elements that have not been included in the 
dataset during the optimization Ell- 

Most important, uncertainties in the EDF could also, 
in principle, be propagated to cases where the EDF is only 
one of several theoretical components, each with a few 
sources of uncertainties. The calculation of low-lying ex¬ 
cited states within the quasiparticle random phase approx¬ 
imation (QRPA) is a straightforward example: it typically 
contains approximation of its own (symmetry restrictions, 
limited model space, etc.), but it is also strongly depen¬ 
dent on the EDF. 

Let us firmly reassert here that in both cases, propa¬ 
gating uncertainties estimated using covariance of Baye¬ 
sian techniques provides information only about the im¬ 
pact of said uncertainties. The procedure does little to 
provide ways to reduce them. In EDF optimization, nu¬ 
merical errors due to basis or mesh truncation can easily 
(at least in principle) be remedied. Statistical and a for¬ 
tiori systematic uncertainties are much more difficult to 
address without a detailed understanding of the nuclear 
many-body problem. 

Most uncertainty propagation reported in the litera¬ 
ture was performed with covariance techniques. This sit- 



Figure 5. (Color online) Comparison between the fission bar¬ 
rier predictions for 240 Pu made with UNEDF1 (solid line), with 
a refit of UNEDF1 including 17 more masses in neutron-rich 
nuclei (dashed line), together with the 90% confidence interval 
(shaded gray area) obtained from Bayesian analysis; from [97l . 

Bayesian techniques have been introduced only recently 
in nuclear theory in general, and EDF optimization in par¬ 
ticular. As a result, in only a couple of cases have these 
methods been applied to the propagation of uncertain¬ 
ties. In 1 the backward-forward Monte-Carlo algorithm 
[124] , which is a particular implementation of Bayesian in¬ 
ference, was used to estimate the statistical uncertainties 
in Skyrme mass models. In [761197] . the full posterior dis¬ 
tribution of the UNEDF1 Skyrme EDF was determined 
in a statistical setting by using Bayesian inference, with 
uniform prior for x and a Gaussian process to emulate the 
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response of the model r](x). The posterior distribution was 
then sampled and used to estimate uncertainties on the fis¬ 
sion barrier of 240 Pu and the position of the two-neutron 
dripline. The large uncertainties on fission barriers visi¬ 
ble in fig. [5] emphasizes the lack of constraints on model 
parameters, which could be caused by an inappropriate 
choice of experimental data and/or too limited a model 
(here Skyrme EDF). 

In addition to the applications mentioned in the previ¬ 
ous section, a few attempts have been made to propagate 
statistical uncertainties from the nuclear EDF to the cal¬ 
culation of observables that involved another model. For 
example, quantifying the impact of neutron skins on the 
electric dipole polarizability or on the weak-charge form 
factor requires calculating the electric dipole response func¬ 
tion, that is, RPA calculations |125H126ill27| . 


5 Conclusions 


Over the past decade, nuclear density functional theory 
has positioned itself as a candidate for a global, compre¬ 
hensive, accurate, and predictive theory of nuclear struc¬ 
ture and reactions. Thanks to the (very recent) introduc¬ 
tion in this field of standard statistical tools such as covari¬ 
ance techniques or Bayesian inference, the statistical un¬ 
certainties associated with the most common energy func¬ 
tionals such as the Skyrme, Gogny, or relativistic EDF 
have been computed rigorously. The propagation of these 
uncertainties to model prediction in nuclei far from sta¬ 
bility has often highlighted the need to substantially im¬ 
prove the constraints on the parameters of the nuclear 
EDF, irrespective of the origin of the functional itself. 
Progress is thus needed in two complementary directions. 
Better rooting of the nuclear EDF in the theory of nu¬ 
clear forces will provide much needed constraints on the 
expected predictive power of the theory. And, this effort 
should go hand in hand with the generalization of sta¬ 
tistical techniques to the problem of EDF optimization, 
and the always-indispensable conversations with the ex¬ 
perimental nuclear physics and data communities. 

On a practical level, an exciting avenue of research 
would be to extend the use of statistical techniques to 
complex problems where the nuclear EDF is one of sev¬ 
eral theoretical tools used. For example, properties of the 
neutron spectrum in neutron-induced fission are currently 
described within the Hauser-Feshbach approach to nu¬ 
clear reactions. Such calculations require fission fragment 
yields, total kinetic energies, and excitation energies of 
the fragments. These quantities, in turn, are currently ob¬ 
tained from either semi-phenomenological models based, 
for example, on Langevin dynamics [12811129] . or from 
fully microscopic time-dependent generator coordinate method 
calculations II3011131111321 . Either way, these dynamical 
calculations depend on the potential energy surface of the 
nucleus in some pre-defined collective space. For the mi¬ 
croscopic approach, this potential energy surface depends 
on which nuclear EDF is used, how the EDF has been 
fitted, and what types of corrections are included [1331IBT1 
I1341I1351I1361I137] . Ultimately, one would therefore wish to 


propagate the uncertainties all the way through this chain 
of “models,” from the nuclear EDF to the fission spec¬ 
trum. 

A related area of future research would be to define 
a comprehensive framework to address uncertainties. In 
this manuscript, we have insisted on the statistical uncer¬ 
tainties, with only a short discussion of numerical errors. 
However, we have also pointed out that all forms of un¬ 
certainties are related to one another: numerical errors 
are not a constant offset in DFT calculations, and thus 
they propagate in a very nonlinear way into the calcu¬ 
lation of the y 2 , which will impact parameter optimiza¬ 
tion and subsequent uncertainty analysis. The particular 
mathematical formulation of the theory (SR-EDF versus 
MR-EDF, HFB approximation only or HFB plus correc¬ 
tions, etc.) also partially determines which observables can 
be reliably computed by the model. For all others, the 
statistical analysis may reveal that some parameters are 
ill-constrained, not because the data is insufficient, but be¬ 
cause the model is not sensitive to it. Moreover, one should 
work toward incorporating experimental uncertainties. In 
the case of the UNEDF2 parameterization, for example, 
both fission isomer excitation energies and single-particle 
states were included in the fit. Yet these quantities are 
model-dependent, and their “experimental” error is rather 
large. In the future, one should try to incorporate this in¬ 
formation in the determination of EDF parameters. 
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